#!/usr/bin/Rscript
# setwd("/home/zhang/BMK/test/cuffdiff")
setwd("F:/cuffdiff/")

read.fpkm <- function(input){
  if(!require(openxlsx)){
    install.packages("openxlsx")
    library(openxlsx)
  }
  
  input <- read.xlsx(as.character(input))
  input <- input[,c(8, 9, 14)]
  
  for(i in 1:(ncol(input) - 1)){
    input[input[,i] == 0, i] <- 0.01
    input[,i] <- log2(input[,i])
  }
  colnames(input)[1:2] <- c("Control", "Treat")
  return(input)
}


# 画箱线图
draw.fpkm.boxplot <- function(input){
  if(!require(ggplot2)){
    install.package("ggplot2")
    library(ggplot2)
  }
  
  stars = ""
  pvalue <- wilcox.test(as.numeric(as.character(input[,1])), as.numeric(as.character(input[,2])))
  pvalue <- pvalue$p.value
  cat(pvalue)
  cat("\n")
  
  if(pvalue < 0.05 && pvalue >= 0.01){
    stars = "*"
  }else if(pvalue < 0.01){
    stars = "**"
  }
  
  input <- cbind(c(input[,1], input[,2]), c(rep("Control", nrow(input)), rep("Treat", nrow(input))))
  colnames(input) <- c("CPM", "Group")
  input <- as.data.frame(input)
  # window specific
  # 改变字体为Times New Roman
  windowsFonts(NRM = windowsFont("Times New Roman"))
  old_theme <- theme_update(axis.text = element_text(family = "NRM", size = 18, color = "black"),
                            legend.key.size = unit(1.2, "cm"),
                            axis.title = element_text(family = "NRM", size = 18),
                            legend.text = element_text(family = "NRM", size = 18))
  # 开始作死
  # 先取出最大值，也就是最高点的y
  y <- max(as.numeric(as.character(input$CPM)))
  # 添加显著线的线本体
  sl <- data.frame(a = c(1,1:2,2), b = c(y + 0.3, y + 0.5, y + 0.5, y + 0.3))
  
  p <- ggplot() + geom_boxplot(data = input, aes(x = Group, y = as.numeric(as.character(CPM)), fill = Group)) + 
    ylab("log2(FPKM)") + geom_line(data = sl, aes(x = a, y = b)) + # 添加显著误差线
    annotate("text", x = 1.5, y = y + 0.6, label = as.character(stars), size = 8)
  return(p)
}


mrna <- read.fpkm("mrna_fpkm.xlsx")
# 所有RNA的cpm的箱线图
tiff("mRNA_fpkm_all_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.fpkm.boxplot(mrna[,1:2])
dev.off()

lnc <- read.fpkm("lnc_fpkm.xlsx")
tiff("lncRNA_fpkm_all_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.fpkm.boxplot(lnc[,1:2])
dev.off()

mrna <- mrna[mrna$significant == "yes", 1:2]
# 差异表达的所有RNA的cpm的箱线图
tiff("mRNA_fpkm_diff_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.fpkm.boxplot(mrna)
dev.off()

lnc <- lnc[lnc$significant == "yes", 1:2]
tiff("lncRNA_fpkm_diff_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.fpkm.boxplot(lnc)
dev.off()

###################################################################################################################
# 单独画microRNA的
read.micro <- function(input, diff = FALSE){
  if(!require(openxlsx)){
    install.packages("openxlsx")
    if(!require(openxlsx)){
      stop("openxlsx install failed")
    }
  }
  
  if(!require(edgeR)){
    source("https://www.bioconductor.org/biocLite.R")
    biocLite("edgeR")
    if(!require(edgeR)){
      stop("edgeR install failed")
    }
  }
  input <- read.xlsx(as.character(input))
  
  if(diff){
    input <- input[abs(input$logFC) >= 1 & input$FDR < 0.05,]
  }
  
  input <- input[,c("control", "treat")]
  
  # 用来将表达量为0的值替换为0.01
  change.zero <- function(input){
    input[input == 0] <- 0.01
    return(input)
  }
  
  input <- cpm(input)
  input <- apply(input, 2, change.zero)
  input <- log2(input)
  return(input)
}

draw.cpm.boxplot <- function(input){
  if(!require(ggplot2)){
    install.package("ggplot2")
    library(ggplot2)
  }
  
  stars = ""
  pvalue <- wilcox.test(as.numeric(as.character(input[,1])), as.numeric(as.character(input[,2])))
  pvalue <- pvalue$p.value
  cat(pvalue)
  cat("\n")
  
  if(pvalue < 0.05 && pvalue >= 0.01){
    stars = "*"
  }else if(pvalue < 0.01){
    stars = "**"
  }
  
  input <- cbind(c(input[,1], input[,2]), c(rep("Control", nrow(input)), rep("Treat", nrow(input))))
  colnames(input) <- c("CPM", "Group")
  input <- as.data.frame(input)
  # window specific
  # 改变字体为Times New Roman
  windowsFonts(NRM = windowsFont("Times New Roman"))
  old_theme <- theme_update(axis.text = element_text(family = "NRM", size = 18, color = "black"),
                            legend.key.size = unit(1.2, "cm"),
                            axis.title = element_text(family = "NRM", size = 18),
                            legend.text = element_text(family = "NRM", size = 18))
  # 开始作死
  # 先取出最大值，也就是最高点的y
  y <- max(as.numeric(as.character(input$CPM)))
  # 添加显著线的线本体
  sl <- data.frame(a = c(1,1:2,2), b = c(y + 0.3, y + 0.5, y + 0.5, y + 0.3))
  
  p <- ggplot() + geom_boxplot(data = input, aes(x = Group, y = as.numeric(as.character(CPM)), fill = Group)) + 
    ylab("log2(CPM)") + geom_line(data = sl, aes(x = a, y = b)) + # 添加显著误差线
    annotate("text", x = 1.5, y = y + 0.6, label = as.character(stars), size = 8)
  return(p)
}


test <- read.micro("microRNA.xlsx")
tiff("microRNA_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.cpm.boxplot(test)
dev.off()

test <- read.micro("microRNA.xlsx", diff = T)
tiff("microRNA_diff_boxplot.tiff", width = 10, height = 7, compression = "lzw", res = 600, units = "in")
draw.cpm.boxplot(test)
dev.off()